Cubic spline interpolation exercise
scientific computing
numerical methods
Large language models (predictably) learn to represent the semantic meaning of sentences.
Cubic spline interpolation exercise local
Imports
Toggle cells below if you want to see what imports are being made.
Data
import json
with open("hr-data.json", "r") as file:
data_raw = json.load(file)
len(data_raw["dataOriginal"])4248
from datetime import datetime
def get_seconds_since_beginning_of_day(date_string: str):
date = datetime.strptime(date_string, "%Y-%m-%dT%H:%M:%S.000Z")
return date.hour * 3600 + date.minute * 60 + date.second
data = [{**item, "time": get_seconds_since_beginning_of_day(item["time"])} for item in data_raw["dataOriginal"]]
data = data[:-1] # drop last point since it overlaps with the next day
data[:5][{'value': 78, 'time': 15},
{'value': 58, 'time': 30},
{'value': 58, 'time': 45},
{'value': 58, 'time': 60},
{'value': 58, 'time': 75}]
Cubic spline interpolation with scipy
# Create a figure
fig = go.Figure()
# Add scatter plot
fig.add_trace(go.Scatter(x=x_train, y=y_train, mode='markers', name='Original points'))
# Add line plot
fig.add_trace(go.Scatter(x=x_test, y=y_test, mode='lines', name='Fitted line'))
# Add titles and labels
fig.update_layout(
title='Scipy interpolation',
xaxis_title='time',
yaxis_title='heart rate'
)
# Show the plot
fig.show()Custom implementation of cubic spline interpolation
def tridiagonal_linear_system_solver(d_lower, d_main, d_upper, b):
"""Solve a tridiagonal linear system given the main, lower, and upper diagonals, as well as the vector b"""
n = len(d_main)
# Forward sweep
for i in range(1, n):
w = d_lower[i - 1] / d_main[i - 1]
d_main[i] -= w * d_upper[i - 1]
b[i] -= w * b[i - 1]
# Back substitution
x = np.zeros(n)
x[-1] = b[-1] / d_main[-1]
for i in range(n - 2, -1, -1):
x[i] = (b[i] - d_upper[i] * x[i + 1]) / d_main[i]
return x
class CubicSplineCustom:
def __init__(self, x, y):
if not np.all(np.diff(x) > 0):
raise ValueError("x values must be in ascending order.")
self.x = np.array(x)
self.y = np.array(y)
self._find_coeffs()
def _find_coeffs(self):
# Find the relevant diagonals and from self.x and self.y assuming natural boundary conditions
d_lower, d_main, d_upper, b, h = self._compute_diagonals_and_b()
# Compute the M_i's for i = 1, n-1 (since M_0 and M_n are assumed to be 0). This will require solving a
# tridiagonal linear system
ms = tridiagonal_linear_system_solver(d_lower, d_main, d_upper, b)
# Compute the coefficients of the cubic polynomials
self.c = self._compute_coefficients_from_second_derivatives(ms, h)
def _compute_diagonals_and_b(self):
x, y = self.x, self.y
h = np.diff(x)
# Compute the diagonals for the tridiagonal matrix
d_lower = h.copy()
d_lower[-1] = 0 # one of the naturals BCs
d_upper = h.copy()
d_upper[0] = 0
d_main = 2 * np.ones(len(x))
d_main[1:-1] *= h[:-1] + h[1:]
b = np.zeros(len(x))
y_diff = np.diff(y)
b[1:-1] = 6 * (y_diff[1:] / h[1:] - y_diff[:-1] / h[:-1])
return d_lower, d_main, d_upper, b, h
def _compute_coefficients_from_second_derivatives(self, ms, h):
coeffs = np.zeros((4, len(self.x) - 1))
coeffs[0, :] = (ms[1:] - ms[:-1]) / (6 * h)
coeffs[1, :] = ms[:-1] / 2
coeffs[2, :] = (self.y[1:] - self.y[:-1]) / h - (ms[1:] + 2 * ms[:-1]) * h / 6
coeffs[3, :] = self.y[:-1]
return coeffs
def _get_index(self, x):
"""Performs binary search"""
low, high = 0, len(self.x) - 1
while low < high:
mid = (low + high) // 2
if self.x[mid] <= x:
low = mid + 1
else:
high = mid
return min(max(low - 1, 0), len(self.x) - 2)
def interpolate(self, x: float):
# Find which polynomial is appropriate and evaluate at x
idx = self._get_index(x)
dx = x - self.x[idx]
return (self.c[0, idx] * dx + self.c[1, idx]) * dx**2 + self.c[2, idx] * dx + self.c[3, idx]
# Example usage:
custom_spline = CubicSplineCustom(x_train, y_train)
print(custom_spline.c.shape)
print(custom_spline.interpolate(5.5)) # Evaluate at x = 5.5(4, 4246)
92.69934779516517
assert cs.c.shape == custom_spline.c.shape
assert np.allclose(cs.c, custom_spline.c)
np.sqrt(((cs.c - custom_spline.c) ** 2).mean())1.1509499204046055e-17
# Create a figure
fig = go.Figure()
# Add scatter plot
fig.add_trace(go.Scatter(x=x_train, y=y_train, mode="markers", name="Original data"))
# Add line plot
fig.add_trace(go.Scatter(x=x_test, y=y_test, mode="lines", name="Fitted line"))
y_test_custom = [custom_spline.interpolate(x) for x in x_test]
assert np.allclose(y_test, y_test_custom)
fig.add_trace(
go.Scatter(x=x_test, y=y_test_custom, mode="lines", name="Fitted line (custom)", line=dict(dash="longdash"))
)
# Add titles and labels
fig.update_layout(title="Scipy and custom interpolations", xaxis_title="time", yaxis_title="heart rate")
# Show the plot
fig.show()